I recently learned about an interesting way of reducing the dimension of a data set. I’d like to illustrate it with a data set of about 550 historical battles in the UK, extracted from Wikidata using a SPARQL query. No doubt some readers will have seen this idea before, but it’s new to me and I found it pretty.

The following plot shows the raw data, with a crossed swords symbol for the location of each battle. It’s possible to see some patterns on the plot (for example, it probably wasn’t much fun to live on the borders of mediaeval Wales) but it could be clearer.

One way to visualise battles over time is to make a 2D graph with time on one axis and location on another axis. But how do you plot a 2D location on a 1D line? One answer is to use a space-filling curve, which in this case will be a curve which snakes all over the map, visiting every location once. This won’t necessarily preserve the relationship between points (points which are close together in 2D could be far apart in 1D) but at least points which are close together in 1D are guaranteed to be close together in 2D.

One example of a space-filling curve is the Hilbert curve. A typical Hilbert curve can be created by writing a function which takes an integer n as input and produces a $2^n \times 2^n$ matrix like this:

This matrix should be interpreted as a curve which starts at zero and follows the numbers in order, like this:

As n increases, the curve becomes more and more complicated. The visualisation below uses a $64 \times 64$ grid and a curve like this:

For the final visualisation, I used a rectangle between $-6$ and $1.8$ degrees longitude and $49.8$ and $58$ latitude. I added the locations of some cities to the plot and plotted squares which were mostly sea as blue lines.

I quite like the visualisation. Firstly, you can see some of the geography. There are three clusters of points which, going from left to right, seem to correspond roughly to Wales, Scotland, and England (the original data contained part of Ireland, but I cut out most of it using the bounding rectangle). Periods of extended conflict over a particular place appear as vertical columns of points. When there was a general conflict over a wide area, you see a horizontal line of points. Some features of the data are much more obvious on this visualisation than on the original map.

I’ve highlighted a few interesting features. The most obvious is the line of points highlighted in yellow, which correspond to the English Civil War. The points highlighted in purple are the Jacobite Rebellion, which up to now were the last major land battles in the UK. The long Conquest of Wales which was also visible on the original map is highlighted in turquoise.

Various other patterns can be detected by staring at the map. The weird-looking row of points near London highlighted in green are various battles from the twentieth century, mostly associated with the two world wars, including the Battle of Britain and various naval engagements.

Overall, I thought this was a very interesting technique, and I wonder where it comes from. Who came up with the idea of using a space-filling curve for dimension reduction? Is there an optimal way to do it (i.e. given some data, is there a “best” way of choosing a space-filling curve to minimise distortion in your plot?) I don’t know!

My recent Dogs of New York, post got circulated a little. On Friday I got a very nice note from William Kornblum, an urban and community studies sociologist who these days is emeritus at CUNY. In the 1980s, Bill studied the people who use Central Park and wrote a report for the Central Park Conservancy. Some of that work counted dog-walkers.

He reminded me that parks in New York City see a lot of dogs. NYC Parks estimates about 1.5 million canine visits to Central Park each year. That is about 4,100 visits a day, to one park. He also asked:

Is there a size gradient perhaps? Larger dogs out on the periphery where there is more room? Smaller dogs in the dense apartment neighborhoods?

This is a good question; one the data I have can almost answer. The main problem is that the license data contained in nycdogs do not record the size or weight of any particular dog. Licenses do record the breed that the owner reports (though, as we saw before the modal breed reported is “Unknown”). But let’s see how far we can get by leveraging our knowledge of breeds.

To begin we need to give each breed a size. For this I used the American Kennel Club’s breed weights. The AKC lists a weight range for each breed it recognizes, from its breed standards. Tiago Filho collected these into one table for 277 breeds in 2020. I had to make a few adjustments to reconcile it with the breeds recorded in nycdogs, most notably for two popular city breeds. First, the AKC lists one Dachshund range that covers both sizes within this breed. Owners often record “Miniature” in the breed name, so I gave those dogs the AKC miniature limit of 11 pounds. Second, Pit Bulls are among the most common dogs in the city, but to my knowledge the AKC does not register them. I used the United Kennel Club’s American Pit Bull Terrier standard instead.

I assigned each breed a “typical weight”, which is the middle of its range. Some standards give only one figure, such as “not more than 6 pounds” for a Chihuahua. In those cases I used that number. This is not ideal. A breed standard describes a show dog, not the dog you see in the park. (Though, of course, in some NYC parks those might be the same. But not in general.) Real dogs will be heavier or lighter than their standard. Still, a Great Dane is reliably bigger than a Yorkie, and it’s the relative differences we’re interested in.

The second problem is that many dogs have no breed we can use. I work with an approximation of one record per dog, as in the first post, which gives about 358,000 dogs. I can only assign a breed size to about 60 percent of them:

• About 22 percent are crossbreeds, such as “Labrador Retriever Crossbreed”. A mix could be any size. I chose not to guess.
• About 10 percent have no breed recorded, mostly “Unknown”.
• About 6 percent are “designer” breeds, such as the up-and-coming Goldendoodles, Maltipoos, and Morkies. No kennel club sets a standard for these.
• About 1.5 percent are Poodles with no size given. A Toy Poodle and a Standard Poodle differ by about 50 pounds, so again, I chose to drop these.

Now, if the missing dogs were spread evenly across the city, the missingness would not matter all that much. We can try to check that. Here is a map showing the share of dogs in each ZIP code that have a breed size.

Coverage runs from about half to about two-thirds of dogs in each ZIP code. It is highest on the West Side of Manhattan and in parts of eastern Queens. It is lowest in central Brooklyn. So, dogs are not MCAR in the city. The maps and analysis below is restricted to describing dogs with a known breed, and that group is not quite the same in every part of the city.

Where the big dogs are

For each ZIP code with at least 50 dogs of known size, I calculated the median typical weight. (As before, I use the city’s MODZCTAs, which are ZIP code areas with some small ones merged.)

I didn’t try to make the weight categories too fine-grained. Again, the weights are derived from the breed standard, not individual dogs, and so the distribution is quite lumpy. Each dog is counted as having the typical weight of its breed. Very popular breeds therefore pile up on single values. Many ZIP codes have a median of exactly 12.5 pounds, for example. (That’s the Shih Tzu’s weight.)

Even so, we still see a clear pattern. The ZIP codes with the biggest dogs are in lower Manhattan and brownstone Brooklyn. Park Slope is at the top, with a median of 29 pounds. Woodlawn and City Island in the Bronx, and the Rockaways in Queens, also have big dogs. The smallest dogs are in most of the Bronx, southeast Queens, and East New York and Brownsville in Brooklyn. There, the median dog weighs 12.5 pounds.

So Bill’s idea is partly supported by the evidence. Some of the city’s edges, where people have more room on average, do have bigger dogs. But, on the other hand, some of the densest apartment neighborhoods in the city have the biggest dogs of all.

Size classes

A median can hide a lot. I sorted the breeds into five size classes. I used the weight bands that you see on bags of dog food: Toy (up to about 9 pounds), Small (9 to 22), Medium (22 to 55), Large (55 to 99), and Giant (over 99). This is kind of arbitrary, and not even arbitrary in the “AKC Standard” sense of arbitrary. But I don’t think it’s unreasonable.

Across the city, here is how the dogs with a known size divide up:

• About a quarter are Toy dogs.
• About a third are Small.
• About a quarter are Medium.
• 14 percent are Large.
• Fewer than 2 percent are Giant.

Here is the share of each class within each ZIP code.

Toy dogs are most common in the Bronx and in parts of Brooklyn and Queens. Medium dogs cluster in somewhat fancier parts of Brooklyn and lower Manhattan. Large dogs are most common on Staten Island, in the Rockaways, and (again) in parts of Manhattan. Giant dogs are rare everywhere.

The same numbers aggregated up to the level of boroughs are a little easier to compare:

Small dogs are about a third of all licensed dogs in every borough. The differences are at the ends. A third of the Bronx’s dogs are Toy dogs, compared with about a fifth of Manhattan’s. Only 9 percent of the Bronx’s dogs are Large, compared with 16 percent in Manhattan and 17 percent on Staten Island. Staten Island, the least densely-populated borough, has the largest share of big dogs.

Mo Money, Mo Doggo

Manhattan and brownstone Brooklyn are dense, but they are also rich. The Bronx is dense too, but much poorer. Happily, thanks to both NYC OpenData and the Census Bureau, we have relevant housing and income measures available that we can construct for each ZIP code.

• Floor area per home. This comes from PLUTO, the city’s database of every tax lot. I added up the residential floor area and the number of homes on every lot in each ZIP code, and divided. This is gross floor area, so it includes halls and stairs. It’s a rough measure of apartment size. You can use my nycpluto package to access these data in R.
• Percent of homes in one- or two-family houses. This is also from PLUTO, using each lot’s building class.
• Median rooms per home, percent of homes that are owner-occupied, and median household income. These come from the Census Bureau’s American Community Survey, using the five-year estimates for 2020 to 2024. You can query the Census API directly for these, or use my nycdemog package to get them.

The Census reports these figures for ZCTAs, its version of ZIP code areas. Most MODZCTAs contain only one ZCTA. Where a MODZCTA contains more than one, I added up the counts. For the medians, I took an average of the ZCTA medians, weighted by the number of households. The Census caps median income at $250,000. A few of the richest ZIP codes in Manhattan hit that cap.

I fit a regression with this data where the outcome is the typical weight of the dogs in each ZIP code. I use an average rather than the median here, because the median is so lumpy. Specifically, I average the logged weights, which is the same as taking the geometric mean. Breed weights are very skewed. The middle dog weighs about 17 pounds, but a few breeds weigh well over 100. And differences in weight are proportional. Five pounds is a big difference between two toy dogs, but not between two retrievers. On the log scale, effects become percent differences in weight. The predictors are the five measures above, plus the borough. There are 177 ZIP codes in the model.

The five measures use different units: dollars, square feet, rooms, and percentages. This makes the coefficients hard to compare directly. What we can do is ask “If we take a ZIP code that is low on this measure (at the 25th percentile of ZIP codes) and move it to a high value (the 75th percentile), holding everything else constant, by what percent do its dogs get bigger or smaller?” For the boroughs, the figure shows the difference between each borough and Manhattan.

Income is the clearest result. A ZIP code with a median household income of $116,000 has dogs that are about 17 percent heavier, on average, than one with a median income of $67,000. The housing measures make much less difference once we know income, and they do not point the way Bill expected. The share of homes in one- or two-family houses and the share of owners have effects close to zero. Median rooms has a small positive effect, maybe, with an interval that includes zero. Floor area per home has a small but clearly negative effect. Moving from about 910 to about 1,100 square feet per home goes with dogs about 3 percent smaller, at the same income. Together, the four housing measures do improve the model a little, by more than we would expect by chance. But mostly that is the floor area effect, and it runs the wrong way for the idea that more room means bigger dogs.

Among the boroughs, Staten Island’s dogs are about 10 percent heavier than Manhattan’s at the same levels of income and housing. That fits the idea of room on the periphery, although the interval is wide, running from about 1 to 20 percent. The other boroughs do not differ much from Manhattan.

In terms of how much variation we can explain, income alone accounts for about 53 percent of the variation in typical dog size across ZIP codes. Adding the borough raises that to 54 percent. Adding all four housing measures raises it to 57 percent.

Here is the income effect on its own. This is a partial regression plot. For both dog size and income, we remove the part that the borough and the housing measures predict. What remains shows the relationship between income and dog size among ZIP codes with similar housing in the same borough.

The relationship is pretty good. Each doubling of median household income goes with dogs about 22 percent heavier. The points come from every borough. So there is a size gradient, as Bill speculated, but it’s not really geography at work here, it’s money. Some of the biggest dogs in New York live in dense, expensive neighborhoods such as Park Slope, SoHo, and the West Village. The smallest dogs live in the Bronx, where incomes are much lower. Once we know a ZIP code’s income, its housing tells us only a little more about its dogs, and not in the expected direction. If anything, more floor area per home goes with slightly smaller dogs. Staten Island is the one place where extra room may add to dog size, beyond what income predicts.

Again, as a paid-up social science rather than infographic guy, I feel like I have to close with the caveats. Remember these are breed weights, not dog weights. Four in ten dogs have no breed size. And, like we said before, these are the licensed dogs only. Licensing rates probably vary with income too. Finally, this is a comparison of ZIP codes, not of households. It tells us that richer ZIP codes have bigger dogs on average. It does not tell us that richer people choose bigger dogs, although that seems likely.

Why would money matter more than space? Big dogs cost more to feed and to treat at the vet. They are harder to keep in buildings that limit pets by weight. And their owners may need to pay someone to walk them while they are at work. Perhaps those many daily dog visits to Central Park are part of the answer too. If you live near a very large park, maybe it’s easier to have a very nice apartment that’s comparatively small for your very large dog.

It’s really easy to make your own agent. You should do it once, just to understand how this technology works.

LLM agents are becoming pervasive, with Meta’s Muse perhaps the first consumer grade agent. They are only going to get more common. If you know a bit of R you can program your own, and it makes the whole thing feel a lot less like magic.

I’ve updated my live book with a new section on making your own agent with R.

The tutorial runs you through using the ellmer package to connect to an LLM. Then you create tools. These are functions that let the agent take actions, like writing and running R code. You tell the agent how to use the functions, then write a parser that checks the LLM response for tool use and executes the tools if it finds any. Then you put all that in a while loop so the agent re-consumes the output of the code it wrote.

That’s it.

This all runs more smoothly if you ‘register’ tools formally, because the agent is less likely to make syntax mistakes when using them. ellmer supports this.

Packages

#install.packages("remotes")
#install.packages("glmnet")
#install.packages("randomForest")

remotes::install_github("thierrymoudiki/caretS3")

caretS3::caretS3 Tabular Machine Learning

 data(mtcars)

 # Prepare data
 X_reg <- mtcars[, -1]  # All except mpg
 y_reg <- mtcars$mpg     # Target variable

 # Split into train/test
 set.seed(123)
 idx_reg <- sample(nrow(X_reg), 0.7 * nrow(X_reg))
 X_train <- X_reg[idx_reg, ]
 y_train <- y_reg[idx_reg]
 X_test <- X_reg[-idx_reg, ]
 y_test <- y_reg[-idx_reg]

 mod <- caretS3::caretS3(X_train, y_train, method = "glmnet")

Registered S3 method overwritten by 'future':
  method               from      
  all.equal.connection parallelly

          Mazda RX4       Mazda RX4 Wag      Hornet 4 Drive             Valiant 
           22.88234            22.41295            21.12871            20.55533 
         Merc 450SE          Merc 450SL Lincoln Continental       Toyota Corona 
           14.37278            15.58185            10.76004            24.86612 
         Camaro Z28    Pontiac Firebird 
           13.55088            15.78322           Mazda RX4       Mazda RX4 Wag      Hornet 4 Drive             Valiant 
           20.09445            20.05342            19.60069            19.34515 
         Merc 450SE          Merc 450SL Lincoln Continental       Toyota Corona 
           15.35807            15.32211            12.98877            24.25831 
         Camaro Z28    Pontiac Firebird 
           15.33895            16.01282 

ahead::dynrmf time series forecasting

par(mfrow=c(2, 2))
plot(ahead::dynrmf(AirPassengers, h=20, fit_func = caretS3::caretS3, fit_params = list(method="glmnet", alpha=0, lambda=0.1)))
plot(ahead::dynrmf(AirPassengers, h=20, fit_func = caretS3::caretS3, fit_params = list(method="glmnet", alpha=0.5, lambda=0.1)))
plot(ahead::dynrmf(AirPassengers, h=20, fit_func = caretS3::caretS3, fit_params = list(method="glmnet", alpha=1, lambda=0.1)))
plot(ahead::dynrmf(AirPassengers, h=20, fit_func = caretS3::caretS3, fit_params = list(method="glm")))

This blog post is about me trying to make a better version of an existing diagram that I think is important, but difficult to read in the version it keeps popping up on my Bluesky feed.

Robin J Brooks of the Brookings Institution has been using a useful diagram, I think since March 2026, to present a set of illustrative theoretical impacts of the USA-Iran war on crude oil prices. For example, he reproduces it in this recent blog post. Relative to many other voices, Brooks has generally been bearish on crude oil prices in this conflict, with his approach very much based on this diagram. The diagram is attributed to Haver Analytics, although I can’t find the original.

I think it’s an important graphic because it does make us confront that the demand curve for crude oil isn’t vertical. There are in fact estimates of the price elasticity of demand for crude oil, and they vary, but certainly sometimes they come out somewhere between -0.1 and -0.2. Usually the price elasticity of demand is introduced to beginner economists as a value of -0.1 means “a price increase of 1% will lead to a decrease in the quantity demanded of 0.1%”. But for this exercise, we will turn the relationship around—if the equilibrium point where supply meets demand is forced down by 0.1% (because of a shift in the supply curve), this implies a price increase of approximately 1%.

In the case of the current war, at its peak, events removed an amount of oil from the market of somewhere around 10% of normal global oil supply. Up to 20% of world supply has been exposed to risk, although I don’t think we ever reached a reduction of that magnitude. The amount removed is much less right now at the time of me writing (meaning, oil is flowing at a much higher rate than a few months ago), but who knows where it will be in a few weeks. The above reasoning implies an impact of a 10% cut in supply on the price of 100% (ie doubling) if the price elasticity of demand is -0.1; or around 50% if it is -0.2. This is basically what the diagram used by Brooks shows.

In doing this, the graphic uses a linear elasticity approximation. I’m not sure this really applies for the scale of change we’re talking about here, with a 10% or 20% drop in supply. But it’s a starting point.

Also note that none of this factors in a change in the demand curve. These are changes up and down the existing demand curve. Over a longer period, the curve will shift with changes in income, policy, buyer attitudes and technology, in particular through the obvious mechanism of accelerated transition to clean, renewable electricity. For now, we are putting those changes aside.

We’re also putting aside how drawing down on existing inventory can (and certainly did) cushion markets from the reduction in supply. We’re trying to illustrate a point and get a plausible range given various assumptions; not actually forecast prices.

The perspective plot

Unfortunately, I find the chart from Haver Analytics very difficult to follow. It’s trying to show the relationship between two key variables—the oil successfully exported from the Gulf, and resulting impact on prices—and one uncertain parameter, the short-term price elasticity of demand for crude oil. We don’t know this elasticity for sure, but it’s plausibly somewhere between -0.1 and -0.2. So we want to see a range of impacts on price given this unknown parameter.

Note: in the plots below, I follow Haver Analytics in showing the absolute magnitude of the price elasticity of demand, dropping the negative sign.

To deal with this, the original uses a 3-dimensional perspective plot. Here’s my rough reproduction in R of the original chart:

I think my version is actually a bit easier to read than Haver Analytics’ original because of the colour shading, but it’s only a marginal improvement. The issue is the common one with these sorts of 3-dimensional plots—it’s just very hard to read across from the flexing surface to the values on the three axes, and it’s even harder to take in a general trend or tendency in the way we find straightforward with a 2-dimensional scatter plot or line plot.

I simulated the illustrative data behind the diagram and made that 3-dimensional perspective plot with this code:

library(tidyverse)
library(plot3D)
library(scales)
library(glue)
library(RColorBrewer)

#---------------Simulating data---------------------------

global_oil_consumption  <- 100
original_gulf_exports <- 20
prewar_oil_price <- 72.5

elasticities <- seq(from = 0.1, to = 0.2, length.out = 11)
exports <- 2:original_gulf_exports
length(exports)
length(elasticities)

d <- expand_grid(
    elasticity = elasticities,
    export = exports
) |>
    mutate(
        lost_oil = (original_gulf_exports - export) / global_oil_consumption ,
        price_rise = lost_oil / elasticity * 100
    )

#---------------common caption and title across all plots-------------------

the_caption <- str_wrap(
    glue(
        "Illustrative calculation assuming global oil consumption of ",
        "{global_oil_consumption } million barrels per day and pre-war Gulf exports of ",
        "{original_gulf_exports} million barrels per day. ",
        "Re-worked by freerangestats.info based on original analysis by Haver Analytics."
    ),
    115
)
the_title <- "Illustrative impact of Gulf export cuts on oil prices"
#----------------Perspective plot---------------
# this recreates the original plot by Haver Analytics

prices_matrix <- d |>
    arrange(export, elasticity) |>
    select(elasticity, export, price_rise) |>
    pivot_wider(names_from = elasticity, values_from = price_rise) |>
    select(-export) |>
    as.matrix()

reds <- colorRampPalette(c("mistyrose", "darkred"))(100)

persp3D(
    x = exports,
    y = elasticities,
    z = prices_matrix,
    col = reds,
    xlab = "\nOil exported from the Gulf",
    ylab = "\n\nPrice elasticity of demand",
    zlab = "\nRise in oil prices (%)",
    main = the_title,
    # box type must be one of “b”, “b2”, “f”, “g”, “bl”, “bl2”, “u”, “n”
    bty = "g",
    ticktype = "detailed",
    colkey = FALSE
)

mtext(
    the_caption,
    side = 1,
    line = 4,
    cex = 0.8,
    col = "grey30"
)

Improvement 1: heatmap

Often, a perspective plot is improved by turning it into some kind of heatmap. These use the same type of data—a grid of values on x and y axes, with a variable of interest on a third dimension. In the perspective plot, that third variable is mapped to the z vertical axis; in a heatmap we map it to the colour of the fill.

Here’s my first attempt at improving the perspective plot by redrawing it as a heatmap:

I think this is already a big improvement, but it’s come at the cost of nearly turning the graphic into a table. This is easier to read, but doesn’t feel like we’re really using the power of visualisation.

Here’s the code that draws the heatmap version:

#-----------------Heat map--------------------
elasticities2 <- c(0.1, 0.12, 0.14, 0.16, 0.18, 0.2)
exports2 <- 1:10 * 2

p1 <- d |>
    filter(round(elasticity, 2) %in% elasticities2) |>
    filter(export %in% exports2) |>
    mutate(text_col = ifelse(price_rise >= 100, "Pale", "Dark")) |>
    ggplot(aes(x = export, y = elasticity, fill = price_rise)) +
    geom_tile() +
    geom_text(
        aes(label = percent(price_rise, scale = 1), colour = text_col),
        size = 4
    ) +
    scale_fill_gradientn(colours = reds, label = percent_format(scale = 1)) +
    scale_y_continuous(breaks = elasticities2) +
    scale_x_continuous(breaks = exports2) +
    scale_colour_manual(
        values = c("Pale" = "white", "Dark" = "black"),
        guide = "none"
    ) +
    labs(
        title = the_title,
        subtitle = str_wrap(
            "Illustrative price increase for crude oil depends on the price elasticity of demand and the amount of oil that is exported from the Gulf.",
            100
        ),
        x = "Oil exported from the Gulf (in million barrels per day)",
        y = "Price elasticity of demand for crude oil",
        fill = "Price increase for crude oil",
        caption = the_caption
    )

print(p1)

Improvement 2: line plot

So I had a go at a line plot. I think the most important relationship here is between the amount of oil successfully exported from the Gulf and the change in world crude oil price. So I make this the x and y axes. The visual relationship between x and y is the most striking thing for many readers of charts; certainly for me. The uncertainty about the price elasticity of demand is of secondary interest, and we can show that just by drawing a different line for a bunch of different plausible values.

Here’s what I came up with:

I think this is definitely the best version.

One of its advantages is that I can use a secondary axis to convert the percentage price increases into an implied dollar value per barrel, based on a pre-war price of around US$72.50 per barrel. This is a nice addition for comprehensibility.

I think the fan of colour-coded lines is also a good way of showing the sensitivity of results to the different possible assumptions about the price elasticity of demand. I’ve deliberately colour-coded these based on the distance from the middle of the range of possible assumptions, rather than from lowest to highest - because I think this better draws attention to this as a type of uncertainty.

But the most important thing with this chart is definitely that we have put the two most important variables on the x and y axes. It means, for example, we can simply see the implied relationship as a trend, which is so hard to see in the original.

It also makes it easy to read across from today’s price and up from the crude exports level (no-one knows what this is today, but say somewhere between 7 and 15 million barrels a day) to hit the fan of different lines, and see a) what any particular value implies about the price elasticity of demand and/or b) if there is an extra price impact being built in, perhaps from other information on inventories, perceived risk, known coming changes in supply and demand, etc.

For what it’s worth, with plausible values for today, the chart implies that there is a bit extra built into the prices. My presumption is that this is likely a (whether rationally founded or not) fear that the exports could well rapidly drop from the current reasonably high levels.

Here’s the code for drawing the line plot version:

#-----------------line plot------------------
annot_col <- "grey12"

p2 <- d |>
    mutate(
        elasticity_distance = abs(elasticity - 0.15),
        is_midpoint = ifelse(round(elasticity, 2) == 0.15, "solid", "dashed")
    ) |>
    ggplot(aes(
        x = export,
        y = price_rise,
        group = elasticity,
        colour = elasticity_distance
    )) +
    geom_line(aes(linetype = is_midpoint)) +
    geom_text(
        data = filter(
            d,
            export == min(export),
            round(elasticity, 2) %in% c(0.1, 0.15, 0.2)
        ),
        aes(label = glue("Elasticity = {elasticity}")),
        hjust = 1,
        position = position_nudge(x = -0.1),
        colour = annot_col,
        size = 2.8
    ) +
    scale_linetype_manual(
        values = c("solid" = 1, "dashed" = 5),
        guide = "none"
    ) +
    scale_y_continuous(
        label = percent_format(scale = 1),
        sec.axis = sec_axis(
            ~ prewar_oil_price * (1 + . / 100),
            name = "Implied crude oil price (US$/barrel)",
            label = dollar
        )
    ) +
    annotate(
        "text",
        x = 17,
        y = 0,
        hjust = 1,
        size = 2.5,
        label = "Pre-war situation",
        colour = annot_col
    ) +
    annotate(
        "segment",
        x = 17,
        xend = 19.5,
        y = 0,
        yend = 0,
        arrow = arrow(length = unit(0.15, "cm")),
        colour = annot_col
    ) +
    expand_limits(x = c(0, 20)) +
    scale_colour_gradientn(
        colours = brewer.pal(9, "Blues")[9:4],
        breaks = c(0, 0.025, 0.05),
        limits = c(0, 0.05)
    ) +
    labs(
        x = "Oil exported from the Gulf (million barrels per day)",
        y = "Price increase for crude oil",
        colour = "Distance from midpoint of assumed elasticity range:",
        title = the_title,
        subtitle = str_wrap(
            "Illustrative price increase for crude oil depends on the price elasticity of demand and the amount of oil that is exported from the Gulf.",
            100
        ),
        caption = the_caption
    ) +
    theme(
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_line(colour = "grey88"),
        panel.grid.minor = element_blank()
    )

print(p2)

That’s all for today. I like my final graphic, but I like the underlying analytical framework even more. I haven’t set out to criticise the original graphic, but to redraw the chart in a way that’s easier for me to read. Turns out the most important thing is to put the two most important variables on the x and y axes.

Do you know the interrupted Goode homolosine projection?

2 min.

What would happen if we peeled the Earth like an orange? The interrupted Goode homolosine projection gives us something close to that! This equal-area projection preserves relative areas while splitting the world map into lobes. Add a few shades of orange and our world is ready to peel.

Let’s make it with sf, giscoR and ggplot2. First, we get the coastlines and draw the outline, following the cuts between the lobes.

library(sf)
library(ggplot2)
library(giscoR)

world <- gisco_get_coastal_lines(resolution = 10)
goode_crs <- "+proj=igh"
lat <- c(90:-90, rep(c(-90:0, 0:-90), 3), -90:90, 90:0, 0:90, 90)
lon <- c(
  rep(180, 181),
  rep(c(80.01, 79.99, -19.99, -20.01, -99.99, -100.01), each = 91),
  rep(-180, 181),
  rep(c(-40.01, -39.99), each = 91),
  180
)
outline <- st_sfc(st_polygon(list(cbind(lon, lat))), crs = 4326) |>
  st_transform(goode_crs)
rectangle <- st_as_sfc(st_bbox(outline) * 1.1)
mask <- st_difference(rectangle, outline)

Next, we draw the meridians and parallels within each lobe so that they do not connect across the cuts.

lobes <- list(
  c(-180, -40, 0, 90), c(-40, 180, 0, 90),
  c(-180, -100, -90, 0), c(-100, -20, -90, 0),
  c(-20, 80, -90, 0), c(80, 180, -90, 0)
)
lines <- lapply(lobes, function(b) {
  west <- b[1] + 0.01
  east <- b[2] - 0.01
  south <- max(b[3], -89.99)
  north <- min(b[4], 89.99)
  meridians <- seq(-180, 180, 30)
  meridians <- meridians[meridians > west & meridians < east]
  parallels <- seq(-60, 60, 30)
  parallels <- parallels[parallels >= south & parallels <= north]
  c(
    lapply(meridians, function(x) {
      st_linestring(cbind(x, seq(south, north, length.out = 181)))
    }),
    lapply(parallels, function(y) {
      st_linestring(cbind(seq(west, east, length.out = 361), y))
    })
  )
})
grat <- st_sfc(do.call(c, lines), crs = 4326) |>
  st_transform(goode_crs)

Finally, we add the orange palette. The mask hides anything outside the outline and keeps the edges of our peel tidy.

ggplot() +
  geom_sf(data = outline, fill = "#FFF0D9", colour = NA) +
  geom_sf(data = world, fill = "#F9B566", colour = "#B34D08", linewidth = 0.25) +
  geom_sf(data = grat, fill = NA, colour = "#E8A45B", linewidth = 0.25) +
  geom_sf(data = mask, fill = "white", colour = NA) +
  geom_sf(data = outline, fill = NA, colour = "#E87516", linewidth = 0.4) +
  coord_sf(
    crs = goode_crs,
    xlim = st_bbox(outline)[c("xmin", "xmax")],
    ylim = st_bbox(outline)[c("ymin", "ymax")],
    expand = FALSE,
    datum = NA
  ) +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "white", colour = NA),
    plot.title = element_text(color = "#E87516", face = "bold")
  ) +
  labs(
    title = "Let's peel the world",
    caption = "Interrupted Goode homolosine projection"
  )

References

The outline and masking approach follows Claus O. Wilke’s Interrupted Goode homolosine example. For more on map projections, see Pyranhia’s Projections cartographiques : le guide pour ne pas y laisser sa santé mentale (in French).

Following this previous post and this discussion on Mastodon, I wrote a function that wraps a combination of functions from the ‘terra’ package to project a SpatRaster or SpatVector layer while automatically centering the projection on that layer, to minimize visual distortion:

projectCentered <- function(x,  # SpatRaster or SpatVector                            proj = "laea"  # +proj parameter for PROJ.4 string                            ) {  # version 1.0, by A. Marcia Barbosa  stopifnot(inherits(x, "SpatRaster") || inherits(x, "SpatVector"))  if (terra::crs(x) == "" || is.na(terra::crs(x))) {    stop("'x' does not have a CRS to project")  }  if (inherits(x, "SpatRaster")) {    a <- terra::as.polygons(x * 0)  }  if (inherits(x, "SpatVector") && length(x) > 1) {    a <- terra::aggregate(x)  } else {    a <- x  }  if (!isTRUE(terra::is.lonlat(a))) {    a <- terra::project(a, "EPSG:4326")  }  centr <- terra::crds(terra::centroids(a))  proj <- tolower(proj)  proj <- sub("+proj=", "", proj, fixed = TRUE)  # in case user input prefix  CRS <- paste(    paste0("+proj=", proj),    paste0("+lat_0=", centr[2]),    paste0("+lon_0=", centr[1])  )  tryCatch(terra::project(x, CRS),           error = function(e)             stop("'proj' is invalid"))}

Here is a use case (compared with the default behaviour), using e.g. Lambert Azimuthal Equal-Area projection:

library(terra)library(geodata)countries <- geodata::world(path = tempdir())map <- subset(countries, countries$NAME_0 %in% c("New Zealand", "Australia"))par(mfrow = c(1, 3))plot(map, main = "Unprojected")plot(project(map, "+proj=laea"), main = "Projected, default-centered")plot(projectCentered(map, "laea"), main = "Projected, map-centered")

This is meant to be used for visualization purposes. For computing e.g. areas or distances, it is better to use unprojected layers, with longitude-latitude coordinates in degrees (making sure that spherical geometry or geodesic algorithms are being used), so that the Earth’s curvature is taken into account.

Note that I just plotted the bare minimum to illustrate what the ‘projectCentered’ function does. Maps can be made to look much nicer — even with just terra + base R plots, with no added package dependencies!

Dear rOpenSci friends, it’s time for our monthly news roundup! You can read this post on our blog. Now let’s dive into the activity at and around rOpenSci!

rOpenSci HQ

Team update: thank you Mark Padgham!

Mark Padgham is leaving rOpenSci at the end of October, but his impact on rOpenSci will remain.

Mark joined rOpenSci as a staff member in December 2019 to help us expand software peer review to statistical topics. However, Mark’s impact on rOpenSci has gone far beyond this one project. Mark has not only shaped how rOpenSci reviews statistical software, but has also built much of the tooling that supports the peer review processes, and R packages which not only support but also strengthen our community. If you have submitted a package to rOpenSci, reviewed or edited one, you have almost certainly used something Mark made!

Read more about his work in our blog post, and consider joining our farewell celebration on Wednesday, October 28th at 15h00 Europe Central.

Champions Program update

Our current 2026-2027 cohort is working on their projects around package development and package peer-review. We celebrated our first cohort meeting (across two different gatherings to accommodate time zones) where participants shared their progress and got feedback from other Champions.

Outreach activities are also under way: Linda Cabrera, assisted by her mentor Luis Verde, conducted her outreach activity by leading a workshop on R Projects with RStudio and GitHub for the Ecuadorian Statistical Society, where she had the opportunity to share information about rOpenSci and our programs with the attendees.

The 2025–2026 cohort wrapped up their journey with a closing community call, Más Allá del Código: muestra abierta de los proyectos de nuestros campeon(a|e)s (Beyond the Code: an open showcase of our Champions’ Projects). At the call, the Champions presented their final projects, from creating, improving, and reviewing R packages to the outreach work that brought those projects to their communities. They also shared what their projects led to, including scholarships, new collaborations, and professional growth. The recording and subtitles, all in Spanish, have been published.

Several Champions and mentors will join LatinR 2026 in Medellín, Colombia, to present their work during the program and participate in HackeaR, the hackathon event where we will work on issues, translations, and Pull Requests on base R and on rOpenSci packages and resources.

If you are around, come and join us!

R-Universe updates

We have made good progress on the migration of Bioconductor to R-Universe infrastructure, see also this presentation by Sean Davis. Among other things, packages from the Bioconductor “experiment”, “workflows” and “books” repositories have been added to bioc.r-universe.dev and bioc-release.r-universe.dev. To support this, R-Universe now has a special category of “data” packages, which get a higher size limit but are only built and checked on Linux.

With about 3000 packages, Bioconductor is by far the largest universe, so along the way we made dozens of performance improvements to the frontend and server, making pages and repositories load faster for every universe. Stay tuned for a dedicated blog post on the state of the Bioconductor transition before the end of the year!

Quinceañera

We continue to uncover the stories of the people behind rOpenSci. This month, Mauro Lepore tells us about his journey with rOpenSci and shares lessons learned that can be applied beyond our community.

We’re also sharing a recap of the synchronous events we held, where we discussed the past, present, and future of rOpenSci with the people who joined us to celebrate these 15 years.

Stay tuned for more rOpenSci’s stories.

Coworking

Read all about coworking!

• Tuesday November 3rd, 09:00 Australia Western (01:00 UTC) “Climate Science in R”, with Steffi LaZerte and co-host Elio Campitelli.
  • Explore how R is used to study the climate
  • Meet co-host, Elio Campitelli, and discuss Climate Science in R
• Tuesday December 8th^*, 14:00 Europe Central (12:00 UTC) “Code Linting in R”, with Steffi LaZerte and co-host Etienne Bacher.
  • Read up on Code Linting and apply some linters to your R code
  • Meet co-host, Etienne Bacher, and discuss code linting in general, or flir and Jarl in particular
    * Note that December coworking is a week later than usual

And remember, you can always cowork independently on work related to R, work on packages that tend to be neglected, or work on what ever you need to get done!

Software

The following three packages recently became a part of our software suite:

• ibger, developed by André Leite together with Marcos Wasiliew, Hugo Vasconcelos, Carlos Amorim, and Diogo Bezerra: Tidyverse-friendly interface to the Brazilian Institute of Geography and Statistics (IBGE) aggregate data API https://servicodados.ibge.gov.br/api/docs/agregados?versao=3. Query aggregates, variables, localities, periods, and metadata from surveys and censuses conducted by IBGE. It has been reviewed by Allan Quadros and Diana García-Cortés.

• rcrisp, developed by Claudiu Forgaci together with Francesco Nattino: Provides tools to automate the morphological delineation of riverside urban areas based on a method introduced in Forgaci (2018) https://doi.org/10.7480/abe.2018.31. Delineation entails the identification of corridor boundaries, segmentation of the corridor, and delineation of the river space using two-dimensional spatial information from street network data and digital elevation data in a projected CRS. The resulting delineation can be used to characterise spatial phenomena that can be related to the river as a central element. It has been reviewed by Mark Padgham and Anatolii Tsyplenkov.

• coevolve, developed by Scott Claessens together with Erik Ringen: Fit Bayesian generalized dynamic phylogenetic models using Stan or JAX. The package allows an arbitrary number of variables with different response distributions to coevolve on a phylogenetic tree via an Ornstein-Uhlenbeck evolutionary process, allowing users to assess directionality between coevolving variables. The optional JAX backend enables high-performance sampling with NumPyro and potential GPU acceleration via reticulate. It has been reviewed by Rob Moss and Krzysztof Bartoszek.

Discover more packages, read more about Software Peer Review.

New versions

The following sixteen packages have had an update since the last newsletter: coevolve (v1.2.1), ruODK (v1.6.3), pkgcheck (v0.3.4), rcrisp (v1.0.0), pkgstats (v0.2.5), openalexR (v3.1.2), webchem (v2.0.0), ReLTER (V3.1.2), rerddap (v1.3.1), ibger (v0.2.1), rgbif (v3.9.0), rglobi (v0.3.5), lakefetch (v0.1.14), visdat (v0.6.1), weatherOz (v3.0.1), and nlrx (v0.4.6).

Software Peer Review

There are seventeen recently closed and active submissions and 4 submissions on hold. Issues are at different stages:

• Three at ‘6/approved’:

  • ibger, Access the IBGE Aggregate Data API from R. Submitted by Andre Leite Wanderley.

  • rcrisp, Automate the Delineation of Urban River Spaces. Submitted by Claudiu Forgaci. (Stats).

  • coevolve, Fit Bayesian Generalized Dynamic Phylogenetic Models using Stan. Submitted by Scott Claessens. (Stats).

• One at ‘5/awaiting-reviewer(s)-response’:

  • RAQSAPI, A Simple Interface to the US EPA Air Quality System Data Mart API. Submitted by mccroweyclinton-EPA.

• Two at ‘4/review(s)-in-awaiting-changes’:

  • camtrapReport, Camera-Trap Report Generator. Submitted by Elham Ebrahimi.

  • nert, Curated Access to TERN Environmental Raster Data. Submitted by Max Moldovan.

• Two at ‘3/reviewer(s)-assigned’:

  • tezr, Access Thesis Metadata from Turkiye’s National Thesis Center. Submitted by Emrah Er.

  • rfastlowess, High-Performance LOWESS Smoothing for R. Submitted by Amir Valizadeh. (Stats).

• Two at ‘2/seeking-reviewer(s)’:

  • corila, Sparse modelling with grouped and correlated features allowing for privileged information. Submitted by Armin Rauschenberger. (Stats).

  • grumpy, Read NumPy .npy and .npz Files. Submitted by Hugo Gruson.

• Seven at ‘1/editor-checks’:

  • ocean3d, Three-Dimensional Marine Spatial Analysis. Submitted by Jay Matsushiba.

  • emodnet.wcs, Access EMODnet Web Coverage Service data through R. Submitted by Maëlle Salmon.

  • qualtdict, Generate Variable Dictionaries and Labelled Survey Data for. Submitted by lyh970817.

  • tarpolyglot, Run Python, Julia, Rust, and C++ Inside targets Pipeline Steps. Submitted by Pierre9344.

  • seasight, Tools for Seasonal Adjustment Diagnostics and Reporting. Submitted by Philipp Wegmueller. (Stats).

  • OptSurvCutR, Optimal Survival Cut-Point Discovery for Time-to-Event Analysis with OptSurvCutR. Submitted by Payton Yau. (Stats).

  • HydraR, Stateful Agentic Orchestration for Scientific Reproducibility. Submitted by Ignatius Pang.

Find out more about Software Peer Review and how to get involved.

On the blog

Software Review

• Celebrating the rOpenSci Community by Mauro Lepore. rOpenSci turns 15 and I’ve been a part of this community for a good chunk of that time. Here I share some memories and the lessons I learned. Other languages: Celebrando a la comunidad de rOpenSci (es).

• Celebrating 15 Years of rOpenSci: Then, Now & Next by Yanina Bellini Saibene, Steffi LaZerte, and Noam Ross. Reflections and ideas shared by participants across the two special gatherings named “rOpenSci Community Celebration: Then, Now & Next”.

• Mark Padgham: Building the Machinery of Software Peer Review by The rOpenSci Team. Mark Padgham is leaving rOpenSci at the end of October, and we want to take this moment to recognize his impact and say thank you for helping to shape rOpenSci.

Calls for contributions

Calls for maintainers

If you’re interested in maintaining any of the R packages below, you might enjoy reading our blog post What Does It Mean to Maintain a Package?.

• suppdata, Grabbing SUPPlementary DATA in R. Issue for volunteering.

• ghql, GraphQL R client. Issue for volunteering.

Calls for contributions

Refer to our help wanted page – before opening a PR, we recommend asking in the issue whether help is still needed.

Package development corner

Some useful information for R package developers.

A Safe Space for Learning How to Make Pull Requests

Jonathan Carroll created and wrote about a package repository that is open to PRs by newbies. It’s a playground for learning how to create pull requests. It’s also an inspiration for package maintainers who want to encourage new contributors!

“My target audience for this repository is people who wish they had somewhere to send a real pull request because they don’t yet feel comfortable participating in open source software.”

Check for global state changes in your testthat tests

Did you know about testthat::set_state_inspector()? If you use it by registering some things you don’t want to change (environment variables for instance), it is run before and after each test, and reports whether there’s any discrepancy in the global state. Here’s how it’s used in vcr:

testthat::set_state_inspector(\() {
 temp_files <- dir(tempdir())
 temp_files <- temp_files[!grepl("^callr", temp_files)]
 temp_files <- temp_files[!grepl("^webfakes", temp_files)]

 list(
 temp_files = temp_files,
 wd_files = dir(),
 vcr_config = the$config
 )
})

Steffi LaZerte noticed the function in Hadley Wickham’s and Jenny Bryan’s workshop materials that Mauro Lepore shared in our Slack workspace.

Using just to run a subset of tests

Speaking of tests… Did you know that testthat’s test function have a filter argument to only run tests that match a given string? Etienne Bacher wrote a nice little post about using the just command-line tool to generate aliases for easily running a subset of tests using this filter.

The post also mentions other use cases of just, e.g. defining just lint to run Jarl on the current directory.

Git hashes

Are you aware that Git plans to change its content hashing algorithm? Here’s a good post by Scott Chacon, the author of the Pro Git book, on why it is important and might be problematic.

Last words

Thanks for reading! If you want to get involved with rOpenSci, check out our Contributing Guide. This guide will help direct you to the right place, whether you want to make code contributions, non-code contributions, or contribute in other ways such as through sharing use cases. You can also support our work through donations.

If you haven’t subscribed to our newsletter yet, you can do so though our signup form. Until it’s time for our next newsletter, you can keep in touch with us through our website, Mastodon, or LinkedIn. See you soon!

When forecasting high-volatility financial assets like Bitcoin (BTC-USD), standard point-prediction models—whether classical ARIMA or deep neural networks—often fall short. A point estimate tells you where the model expects the price to go, but it completely ignores market regime shifts and changing risk levels.

In this article, we construct an end-to-end forecasting pipeline in R that combines Machine Learning (Multi-Layer Perceptron via brulee and tidymodels) for directional trend prediction with Financial Econometrics (GJR-GARCH via rugarch) to generate dynamic, volatility-adjusted confidence bounds.

The Architecture Overview

Instead of attempting to forecast level prices directly—which introduces non-stationarity issues—our pipeline operates on log returns and reconstructs price trajectories dynamically.

Key Architectural Decisions & Methodological Rationale

1. Modeling Log Returns Instead of Level Prices

Directly predicting raw prices often creates spurious regression results because price series are non-stationary. By transforming the price series to log returns:

This transformation achieves stationarity, making the target variable suitable for neural network training. Factual price predictions are reconstructed at test time via cumulative exponentiation:

2. Deep Learning for Trend: Multi-Layer Perceptron (brulee)

We use a 2-layer MLP architecture (32, 16) initialized via the tidymodels framework.

• Lag Features: We supply 20 days of lagged log returns to capture short- to medium-term autocorrelation.
• Normalizations & Dropping Zero Variance: Handled within a clean recipes pipeline to prevent data leakage

3. Econometric Volatility Modeling: GJR-GARCH(1,1) on Neural Network Residuals

Machine learning models capture the conditional mean, but their residuals, $\epsilon_t = r_t – \hat{r}_t$ , still contain volatility clustering—periods of high turbulence followed by high turbulence, and calm periods followed by calm.

To model this conditional variance ($\sigma_t^2$), we fit a GJR-GARCH(1,1) model on the MLP training residuals:

• Why GJR-GARCH over Standard GARCH? Crypto markets exhibit leverage effects: price drops often spark higher panic volatility than equivalent price rallies. The indicator term $I_{t-1}$ captures this asymmetry ($\gamma$).
• Why garchOrder = c(1, 1)? The order (1,1) specifies a single lag for the ARCH component ($\epsilon_{t-1}^2$) and a single lag for the GARCH component ($\sigma_{t-1}^2$). Empirical literature demonstrates that a single lag captures the vast majority of volatility persistence without overfitting.
• Why Skewed Student-t (sstd) Distribution? Bitcoin returns violate Gaussian normal distribution assumptions. They display fat tails (extreme events occur frequently) and skewness (asymmetric distribution). The sstd distribution explicitly models both heavy tails and skewness.
• Solver Selection: We use solver = "hybrid" in rugarch to combine optimization algorithms, ensuring convergence even under complex distribution parameters.

When the market is quiet, these confidence bounds contract automatically. During volatile regimes, the bounds expand, giving a realistic probabilistic envelope rather than static, uninformative confidence intervals.

Complete R Implementation

Below is the complete, self-contained R script implementing the full pipeline using tidymodels, brulee, rugarch, and plotly.

# Integrated Forecasting Pipeline with tidymodels, yardstick, and plotly
library(tidymodels)
library(timetk)
library(brulee)
library(modeltime)
library(tidyquant)
library(plotly)
library(scales)
library(glue)
library(ggtext)
library(rugarch) # Added for GJR-GARCH volatility modeling

# 1. Fetch Data & Calculate Log Returns
df_raw <- tq_get("BTC-USD") %>%
  select(date, close) %>%
  filter(date >= last(date) - months(24)) %>%
  drop_na() %>%
  mutate(log_ret = log(close / lag(close))) %>%
  drop_na()

# 2. Add Lag Features on Log Returns
df_prepared <- df_raw %>%
  tk_augment_lags(log_ret, .lags = 1:20) %>%
  drop_na()

# 3. Time Series Split (10 Assessment Window)
split <- time_series_split(
  df_prepared,
  assess = 10,
  cumulative = TRUE
)

df_train <- training(split)
df_test  <- testing(split)

# 4. Recipe Setup (Predict Log Return, Not Level Price)
recipe_mlp <- recipe(log_ret ~ ., data = df_train) %>%
  step_rm(date, close) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_zv(all_predictors()) %>%
  step_normalize(all_numeric_predictors())

# modeltime::naive_reg() requires the date column as a predictor to parse the time index
recipe_naive <- recipe(log_ret ~ date, data = df_train)

# 5. Model Specifications
model_spec_mlp <- mlp(
  hidden_units = c(32, 16),
  epochs = 100,
  dropout = 0.05,
  learn_rate = 0.005
) %>%
  set_engine("brulee") %>%
  set_mode("regression")

model_spec_naive <- naive_reg() %>%
  set_engine("naive")

# 6. Workflow Construction & Model Fitting
workflow_mlp <- workflow() %>%
  add_recipe(recipe_mlp) %>%
  add_model(model_spec_mlp) %>%
  fit(df_train)

workflow_naive <- workflow() %>%
  add_recipe(recipe_naive) %>%
  add_model(model_spec_naive) %>%
  fit(df_train)

# 7. Predictions & Price Reconstruction (Reference Code Logic)
preds_mlp_ret   <- predict(workflow_mlp, df_test)
preds_naive_ret <- predict(workflow_naive, df_test)

df_results <- df_test %>%
  mutate(
    yesterday_close = lag(close, default = first(close)),
    actual      = close,
    pred_brulee = yesterday_close * exp(preds_mlp_ret$.pred),
    pred_naive  = yesterday_close * exp(preds_naive_ret$.pred)
  )

# ==============================================================================
# 7.5. GJR-GARCH (Skewed Student-t) Dynamic Volatility & Confidence Bands
# ==============================================================================
# Calculate MLP residuals on the training set
train_preds_mlp <- predict(workflow_mlp, df_train)$.pred
train_residuals <- df_train$log_ret - train_preds_mlp

# GJR-GARCH(1,1) specification with Skewed Student-t distribution for fat-tails and asymmetry
garch_spec <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
  distribution.model = "sstd"
)

# Fit GJR-GARCH model on training residuals
garch_fit <- ugarchfit(spec = garch_spec, data = train_residuals, solver = "hybrid")

# Forecast conditional volatility for the 10-day test horizon
garch_fore <- ugarchforecast(garch_fit, n.ahead = 10)
dynamic_sigmas <- sigma(garch_fore) # 10-day dynamic conditional standard deviations

# Attach dynamic standard deviations to test results and compute dynamic confidence bounds
df_results <- df_results %>%
  mutate(
    conditional_sigma = as.vector(dynamic_sigmas),
    # Dynamically scaled confidence bounds in price space (± 1.96 * conditional_sigma)
    conf_hi = pred_brulee * exp(1.96 * conditional_sigma),
    conf_lo = pred_brulee * exp(-1.96 * conditional_sigma)
  )

# 8. Accuracy Table Construction (Yardstick on Original Price Scale)
accuracy_table <- tibble(
  .model_id   = c(1, 2),
  .model_desc = c("BRULEE MLP (32, 16) + GJR-GARCH", "NAIVE"),
  rmse        = c(rmse_vec(df_results$actual, df_results$pred_brulee),
                  rmse_vec(df_results$actual, df_results$pred_naive)),
  rsq         = c(rsq_vec(df_results$actual, df_results$pred_brulee),
                  rsq_vec(df_results$actual, df_results$pred_naive)),
  mape        = c(mape_vec(df_results$actual, df_results$pred_brulee),
                  mape_vec(df_results$actual, df_results$pred_naive)),
  mae         = c(mae_vec(df_results$actual, df_results$pred_brulee),
                  mae_vec(df_results$actual, df_results$pred_naive))
)

# Print Model Comparison Table
cat("=== MODEL ACCURACY EVALUATION TABLE (PRICE SCALE) ===\n")
cat(knitr::kable(accuracy_table, digits = 4, format = "simple"), sep = "\n")

# 9. Modern Interactive Plotly Visualization (Reference Style)
brulee_mape <- accuracy_table$mape[1]
df_eval <- df_results

# Custom Hover Text Data Frames
df_plot_actual <- df_eval %>% select(date, actual)  %>% mutate(text_actual = glue("<b>Actual Price:</b> {scales::dollar(actual, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_pred   <- df_eval %>% select(date, pred_brulee) %>% mutate(text_pred   = glue("<b>AI Prediction Price:</b> {scales::dollar(pred_brulee, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_hi     <- df_eval %>% select(date, conf_hi) %>% mutate(text_hi     = glue("<b>GJR-GARCH Upper Bound:</b> {scales::dollar(conf_hi, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_lo     <- df_eval %>% select(date, conf_lo) %>% mutate(text_lo     = glue("<b>GJR-GARCH Lower Bound:</b> {scales::dollar(conf_lo, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))

# Build Base GGPlot Layer
p <- ggplot() +
  geom_ribbon(data = df_eval, aes(x = date, ymin = conf_lo, ymax = conf_hi), fill = "#808080", alpha = 0.18) +
  geom_point(data = df_plot_hi, aes(x = date, y = conf_hi, text = text_hi), color = "transparent", alpha = 0, size = 3) +
  geom_point(data = df_plot_lo, aes(x = date, y = conf_lo, text = text_lo), color = "transparent", alpha = 0, size = 3) +
  geom_line(data = df_plot_actual, aes(x = date, y = actual), color = "#2c3e50", linewidth = 1.2) +
  geom_point(data = df_plot_actual, aes(x = date, y = actual, text = text_actual), color = "#2c3e50", size = 2) +
  geom_line(data = df_plot_pred, aes(x = date, y = pred_brulee), color = "#e74c3c", linetype = "dashed", linewidth = 1.2) +
  geom_point(data = df_plot_pred, aes(x = date, y = pred_brulee, text = text_pred), color = "#e74c3c", size = 2) +
  scale_y_continuous(labels = dollar_format(accuracy = 1)) +
  labs(
    x = "", y = "",
    title = paste0(
      "<span style='color:#2c3e50'>BTC-USD Prices</span> vs ",
      "<span style='color:#e74c3c'>MLP + GJR-GARCH Dynamic Forecast</span><br>",
      "<span style='font-size:12px; color:#555555;'>10-Day Horizon | Out-of-Sample MAPE: ", round(brulee_mape, 2), "</span>"
    )
  ) +
  theme_minimal() +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    plot.background = element_rect(fill = "#ffffff", color = NA),
    panel.background = element_rect(fill = "#ffffff", color = NA),
    panel.grid.minor = element_blank()
  )

# Typography Settings & Plotly Output
font_family <- list(family = "Roboto Slab, Sans-Serif", size = 16)
label_font  <- list(font = list(family = "Roboto Slab, Sans-Serif", size = 13))

ggplotly(p, tooltip = "text") %>%
  style(hoverlabel = label_font) %>%
  layout(font = font_family)

Conclusion

Combining machine learning with traditional econometrics produces a hybrid system that overcomes the limitations of both fields:

1. Machine Learning captures non-linear relationships and trend movements.
2. Econometrics (GARCH) supplies explicit statistical parameters for uncertainty and volatility clustering.

The resulting output provides realistic, dynamic risk envelopes that evolve with changing market regimes—giving quantitative traders and analysts a robust framework for volatility-aware decision-making.

Mark Padgham is leaving rOpenSci at the end of October, and we want to take this moment to recognize his impact and say thank you for helping to shape rOpenSci.

Mark joined rOpenSci as a staff member in December 2019 to help us expand software peer review to statistical topics. However, Mark’s impact on rOpenSci has gone far beyond this one project. Mark has not only shaped how rOpenSci reviews statistical software, but has also built much of the tooling that supports the peer review processes, and R packages which not only support but also strengthen our community. If you have submitted a package to rOpenSci, reviewed or edited one, you have almost certainly used something Mark made!

Making a Mark

Mark’s most visible contribution is to statistical software peer review. Reviewing a package that implements statistical methods is a different job from reviewing a package that accesses an API. Mark tackled the task of translating what “good” means in a statistical package and how a reviewer can assess it. He recruited a new set of editors with the expertise to tackel this work and turned insights into a concrete set of statistical software standards and the guide that accompanies them (rOpenSci Statistical Software Testing and Peer Review Community Call). Mark also built the machinery to go with this new type review: pkgstats gathers information about a package’s structure, helps authors and editors see what a submission contains and how it is organised; srr allows package authors to document their adherence to statistical software standards.

The result of all this work was our first peer-reviewed statistical packages, published in November 2022, and to date we’ve had over 26 statistical packages go through this review! Anyone who has gone through statistical software review has benefited from the groundwork Mark laid.

The engine of software review

Peer review at rOpenSci is a human process, but it runs on a lot of automation, and a good part of that automation is Mark’s work.

Mark built and maintained pkgcheck, which runs automated checks on submitted packages, and the GitHub Action that allows anyone to run the same checks on their own repository. pkgcheck uses goodpractice, another package maintained by Mark, that recently rolled out a host of updates and extension. Mark maintained Buffy, the bot behind our chat review system on GitHub, which lets editors assign reviewers, track deadlines and move submissions along with simple comments, and many other pieces of infrastructure that make the whole system work.

Review is a human process

In early 2025, Mark became our new Software Review Lead. In this role he was in charge of editorial management, curating the current editorial board and planning its future. Mark helped guide peer-review through new challenges, such as how review should handle AI-assisted contributions, resulting in rOpenSci’s first AI policy; or how to deal with large software submissions, resulting in our first group peer review.

Tools for everyone

Mark’s care for rOpenSci goes beyond software. He researched alternatives to GitHub, sharing the results in Beyond GitHub, a thoughtful look at what it would mean for a community like ours to depend less on a single platform.

He has also served on our Code of Conduct committee, work that helps ensure rOpenSci remains a welcoming and safe community.

Mark also builds tools that help the wider R community. The allcontributors package, inspired by the allcontributors project, makes it easy to recognise everyone who contributes to a project, not just those who write code. His work on repometrics explores how to understand the health and activity of software repositories, and his work on deposits looks at how research data and software can be deposited in trusted repositories.

Thank you, Mark

Mark’s work has shaped how we review software, support contributors, and build a more reliable and sustainable open science ecosystem. Thank you, Mark, for building so much of the infrastructure and process that we rely on, and for leaving it in better shape than you found it.

If you have benefited from Mark’s work, reach out to Mark and let him know! You can keep following his work through his website and social media.

Je dirige seul Qognito, une entreprise de conseil et de formation. Pour réaliser certaines activités de mon cœur de métier, je coordonne plusieurs équipes d’agents IA, mes « squads ». C’est ma manière de faire fonctionner une firme minimale : déléguer une partie du travail tout en restant responsable des orientations, des livrables et des engagements pris auprès de mes clients.

La veille sur l’IA alimente ce travail. Une note peut m’aider à actualiser une formation, à examiner le besoin d’un prospect ou à préparer une proposition. Sa valeur tient à ce que je peux en faire ensuite.

Je conserve ces notes dans Obsidian, une application de gestion de connaissances. Des liens internes, les wikilinks, les relient entre elles : je peux retrouver une analyse et la rapprocher d’autres sujets. Ces liens facilitent l’exploration ; ils ne garantissent pas que plusieurs notes apportent des preuves indépendantes. Une erreur peut ainsi devenir un appui pour un raisonnement ultérieur.

J’avais organisé la production avec soin. Un agent prépare la source, un autre rédige l’analyse, un troisième évalue la note. Je décide de sa publication dans ma base. Voici le circuit attendu :

Ce schéma décrit les responsabilités prévues. Il reste à vérifier que les outils respectent ce circuit.

Dans ma formation SAGA-IA, je demande aux décideurs d’examiner ce qu’ils délèguent, les preuves disponibles et le coût du contrôle. J’ai soumis mon propre dispositif à ces exigences. L’audit a commencé par une découverte très simple.

Les règles pouvaient changer avec celui qui les appliquait

Les fichiers qui définissaient les consignes, les critères du juge et les contrôles étaient modifiables par les agents. Celui qui produisait une note avait donc techniquement la possibilité de changer les critères utilisés pour l’évaluer.

Je n’avais pas constaté qu’un agent l’avait fait. Mais la séparation des rôles reposait sur une consigne que les droits d’accès ne faisaient pas respecter.

J’ai fait protéger ces fichiers. Un agent peut proposer une modification ; sa mise en place exige une intervention humaine avec les droits d’administration. Dix tentatives de modification, par des voies différentes, ont ensuite été refusées. Les règles sont restées intactes.

Ces premières corrections précédaient la campagne consacrée aux protections, appelée « Phase A » dans mon audit. Elles avaient déplacé la question : les règles tenaient mieux, mais que pouvait-il encore arriver aux documents sur lesquels je m’appuyais ?

Publier une note ne suffisait pas à la préserver

Au départ, une note rejoignait la base dès sa rédaction, avant le verdict du juge. J’avais ensuite instauré une étape d’attente, puis une autorisation humaine de publication. Le premier cycle réel avait suivi ce circuit : neuf sources traitées, neuf notes examinées et publiées par mes soins.

Pourtant, une note déjà publiée pouvait encore être modifiée. Une source aussi. Le contrôle existant ne repérait pas ces changements. Une approbation donnée hier ne disait donc pas ce que je retrouverais demain.

La Phase A a traité ce problème. Sur une copie de travail, trois altérations volontaires — une note modifiée, une note ajoutée hors circuit et une transcription modifiée — ont été détectées par un nouveau contrôle comparant les documents à une référence conservée. Des verrous ont également été posés sur les notes publiées et les sources préparées.

J’ai accepté une discipline pour moi aussi : plus de retouche directe dans Obsidian. Les changements doivent passer par le circuit prévu et laisser une trace. Sinon, comment distinguer ma correction d’une modification anormale ?

Préserver le contenu ne suffit toutefois pas à le rendre juste. Lors des premiers essais, le juge avait relevé « usually » traduit par « toujours », alors que les contrôles automatiques n’avaient rien signalé. Une nuance devenue certitude peut changer une décision.

Le programme compare les documents ; le juge examine le sens ; je conserve la décision de publication. Cette répartition rend les responsabilités identifiables. Elle ne dispense pas de vérifier la qualité du jugement.

J’ai dû examiner mes propres examens

Un banc d’essai est une répétition sur des cas choisis, dont certains contiennent volontairement des pièges. J’en avais utilisé pour éprouver la veille. Deux constats m’ont obligé à revoir la confiance accordée à leurs résultats.

Le premier concernait le corrigé : les réponses attendues restaient accessibles aux agents. Imaginez un examen dont les candidats pourraient ouvrir les solutions. Cela fragilise l’évaluation, même sans preuve qu’ils les aient consultées. L’accès au corrigé a été protégé, puis son refus éprouvé sur le poste.

Le second concernait les permissions. La copie utilisée pour certains essais était plus restrictive que l’environnement réel. Une action pouvait y demander un accord alors qu’elle passait librement au quotidien. Les résultats décrivaient bien les essais réalisés ; ils ne suffisaient pas à établir la même protection en usage réel.

La dernière étape de la Phase A devait donc réunir les protections dans une épreuve avec les permissions de production. Il fallait observer les tentatives, identifier ce qui les arrêtait et vérifier leurs effets sur les documents.

Vingt tentatives, et une limite qui reste visible

Le 19 octobre, vingt tentatives ont été exécutées et tracées dans les deux environnements d’agents utilisés. Les essais de modification de documents protégés et d’accès au corrigé ont rencontré des refus.

Deux créations de fichiers hors du circuit de publication ont toutefois réussi dans l’un des environnements. Le contrôle les a détectées ; elles ont ensuite été supprimées. La vérification finale a confirmé le retour de la base à son état initial.

Vingt tentatives documentées ne signifient donc pas vingt actions empêchées. Dans le fonctionnement courant, un ajout hors circuit peut rester présent jusqu’au contrôle déclenché par la prochaine publication. Cette limite reste ouverte.

La Phase A est close, avec une réserve précise : le prochain cycle réel doit encore montrer qu’une publication autorisée aboutit à une note effectivement verrouillée. Les refus constatés pendant les essais ne remplacent pas cette vérification de bout en bout.

L’évaluation de qualité menée en parallèle a donné sept cas réussis sur dix, avec une fausse alerte de fidélité et deux cas dont la portée restait à trancher. Protéger les documents n’établit pas que leur contenu s’est amélioré. Un examen rétrospectif des notes publiées reste prévu.

Le temps que le chronomètre ne voyait pas

Pour comprendre le coût de cette maîtrise, il faut revenir au premier cycle réel. L’outil avait enregistré une à trois secondes par publication. J’avais pourtant consacré neuf minutes à relire les neuf notes avant de l’ouvrir. Ce temps a dû être déclaré séparément.

Ces neuf minutes suffisaient à ma décision de publier ; elles ne constituaient pas une comparaison complète de chaque note avec sa source. Une validation humaine doit préciser ce qui a réellement été vérifié et avec quel temps disponible.

La même exigence vaut pour les corrections. Pendant les premiers travaux, j’ai accepté un défaut de recoupement compté deux fois : une seule occurrence avait été constatée, et le script envisagé pour la traiter aurait été coûteux et lui-même faillible. Chaque anomalie ne justifie pas une règle supplémentaire.

J’ai aussi fixé des critères d’arrêt, des responsables et des conditions de réexamen. Ces dispositions devront vivre dans l’usage. À ce stade, je n’ai démontré aucun gain global de rentabilité : la relecture, les reprises, la maintenance et l’audit font partie du coût à apprécier.

Ce que je peux désormais justifier

Le chemin parcouru tient en trois étapes :

Les protections testées apportent des preuves circonscrites. La dernière étape reste à observer.

Je peux mieux expliquer ce que je confie aux agents, ce qui résiste aux tentatives testées et ce qui demande encore un contrôle. Ma responsabilité de dirigeant consiste aussi à décider quelles incertitudes sont acceptables et lesquelles doivent limiter la délégation.

Si vous coordonnez vous aussi une activité avec des agents, mon livre blanc gratuit La facture fantôme de l’IA propose une première lecture sur la supervision, les coûts complets et l’usage du temps gagné.

Pour travailler ces décisions, SAGA-IA — Volume 1 réunit le diagnostic de délégation et l’évaluation de la viabilité économique, avec des cas et des exercices. C’est le prolongement de ce REX : apprendre à justifier ce que l’on engage, les preuves que l’on demande et les conditions auxquelles on poursuivra.

I often need to get mainland from country / administrative region polygons, so I finally decided to write a function to help automate it. It selects the N largest feature(s) in a polygons or lines map. Here’s the code:

selectLargest <- function(x,  # terra SpatVector polygons or lines                          n = 1,  # number of largest features to retain                          dissolve = FALSE,  # whether to group contiguous features                          singleparts = TRUE  # whether to disaggregate multi-part features) {  stopifnot(inherits(x, "SpatVector"),            terra::geomtype(x) %in% c("polygons", "lines"))  if (dissolve) {    x <- terra::aggregate(x)  }  if (singleparts) {    x <- terra::disagg(x)  }  if (terra::is.polygons(x)) {    size_rank <- order(terra::expanse(x), decreasing = TRUE)  } else if (terra::is.lines(x)) {    size_rank <- order(terra::perim(x), decreasing = TRUE)  }  x[size_rank[1:n], ]}

The dissolving of internal borders is optional. Here are some usage examples:

library(terra)library(geodata)# import a USA map and get the conterminous and continental parts:usa <- geodata::gadm("United States", level = 1, path = tempdir())mainland_usa <- selectLargest(usa, dissolve = TRUE)contin_usa <- selectLargest(usa, dissolve = TRUE, n = 2)  # two largest polygons, to include Alaskaplot(usa, main = "USA")plot(mainland_usa, main = "Conterminous USA")plot(contin_usa, main = "Continental USA")# import a UK map and get Great Britain:uk <- geodata::gadm("United Kingdom", path = tempdir())great_britain <- selectLargest(uk, dissolve = TRUE)plot(uk, main = "UK")plot(great_britain, main = "Great Britain")# import a Portugal + Spain map and get the Iberian Peninsula:ptes <- geodata::gadm(c("Portugal", "Spain"), path = tempdir())iberia <- selectLargest(ptes, dissolve = TRUE)plot(ptes, main = "Portugal and Spain")plot(iberia, main = "Iberian Peninsula")

As usual, feedback is welcome!

TLDR: For the life of me, I couldn’t get Ctrl+Shift+M to work for inserting pipes (|>) when editing a Quarto document in Positron visual mode. After countless hours fiddling with keybindings.json, I developed a functional, but obtuse workaround for Windows. It’s not pretty, but it works. This post outlines the approach in case it helps others facing the same problem.

Note: this problem is already on the radar of the development team, so will no doubt be fixed soon. I’ve also only outlined a workaround for Windows. Feel free to reach out if you’ve found a better solution (or a solution that works on other platforms) and I’ll add this to the post.

The problem

I’m a big fan of Quarto as it provides a seamless way to mix data, code and narrative text without ever leaving my IDE. I’ve also been trying to move from RStudio to Positron as I find it snappier, works better with AI tools and works well with Python. Unfortunately, Positron is still experiencing growing pains as a newish IDE, including some shortcuts not working when editing Quarto documents in visual mode.

This isn’t a new issue (e.g: 1, 2 and 3), so is sure to be fixed in the not-too-distant future (likely just after I post this). But, after countless hours messing around with keybindings.json and hassling people from Posit (apologies, Isabella) I asked Claude if it had a solution.

It did.

And it works.

So, I thought I’d share it here, despite the fact it doesn’t spark joy with the minimalist inside of me:

Duct tape and thumbtacks: intercept the shortcut before it reaches Positron

As far as I can tell, the problem comes from Positron / VSCode intercepting some shortcuts when Quarto is in visual mode, which is probably why adjusting the settings did nothing for me. As a (temporary) solution to this, Claude suggested installing a program called AutoHotkey v2 and creating a script operates when Positron is active.

Obviously, as an independent-minded person I did what I’m told and here’s how you can do the same:

1. Install AutoHotkey v2

Get Version 2 from autohotkey.com (the code below will not run under v1).

2. Create the script

Open Autohotkey2 (or AutoHotKey Dash) and select New Script to create a new blank .ahk script in a secluded folder of your choice. Give it a name like ‘positron-keys.ahk’ so future you knows what it’s doing there and use a text editor to replace the contents of the file with the code below (or the equivalent depending on your shortcut preferences):

(SendText sends the literal characters to the Quarto document, so nothing gets re-interpreted by Positron, while #HotIf restricts the behavior to Positron, so it doesn’t interfere with other applications.)

#Requires AutoHotkey v2.0
#SingleInstance Force

#HotIf WinActive("ahk_exe Positron.exe")

^+m::SendText(" |> ")    ;Ctrl+Shift+M -> native pipe
!-::SendText(" <- ")     ;Alt+-        -> assignment arrow

#HotIf

3. Run it

To run the code and launch AutoHotkey, open the ‘positron-keys.ahk’ file. If it worked, a green H should appear in the system tray.

4. Test it

To test it you can now open Positron and try using Ctrl + Shift + M to insert a pipe in a Quarto when in visual mode. Hopefully it works, as this is the only idea covered in this post.

5. Make it persistent

To avoid needing to load the code each time you want to use Positron, you can then add a shortcut to the script that loads at startup. For Windows, you can open the startup folder by typing Win+R and executing the command ‘shell:startup‘. Right-drag positron-keys.ahk into the folder and choose Create shortcut here.

Note: if you have other solutions, please let me know. I’ll happily update this post if others find it useful.

AI use: Claude suggested this solution and wrote the first draft of this post, which I heavily edited.

The post Getting the pipe shortcut to work in Positron appeared first on Giles.

In June 2026, rOpenSci marked its 15th anniversary with three community celebration events: a themed coworking event and two extra events. This post captures the reflections and ideas shared by participants across the two special gatherings named “rOpenSci Community Celebration: Then, Now & Next” that were held to accommodate different time zones.

How the events worked

Our goal was to create a space for people to have conversations, get to know each other and also reflect on rOpenSci’s past, present and future.

After a welcome and a brief introduction, participants joined three successive breakout rounds of 15 minutes each, one for each theme: the past, the present, and the future. Between rounds, attendees took part in a Scavenger Hunt (a regular feature of our co-working sessions), finding something in their physical surroundings that matched a given prompt.

The Scavenger Hunt is a fun moment together, and we also used it as a natural prompt for introductions at the start of the next breakout room — “what did you find?” becomes a nice and easy opening for conversations.

As in our community calls, workshops and co-working sessions, a shared document served as the place where everyone could take notes in real time, add thoughts, and keep contributing for a couple of days after the event ended. Noam Ross named this “process is product” approach, caring about how things happen, not just what happens, as one of rOpenSci’s defining traits.

Round 1 – The Past: What made rOpenSci special?

The first breakout asked participants how they first connected with rOpenSci and what their favourite memory was.

Several people traced their connection back to a talk about rOpenSci or to a personal introduction from a community member who made them feel welcome from the start.

The unconferences were remembered with particular warmth and excitement: the chance to sit next to people, see how they work, talk about R, and be part of something as it was taking shape. “Meeting my heroes”, as one participant put it.

Other things people named as favourite:

• The peer-review process: praised by editors, reviewers and authors alike, both for the quality of the process and for the connections it created.

Being invited to write a blog post after a first package review … made me feel very cool and a valued part of the community.” – Alex Koiter

• The dev guide: a detailed and evolving documentation with knowledge that isn’t easily found elsewhere.

• The Champions Program and the mini-hackathons, spaces that provide friendly and guided entry points for open source contributors and R package developers.

• Co-working sessions are a place were you can find “your people,” a place where you don’t need to explain what R is.

It’s a satisfying community ‘distraction/procrastination’ from our normal jobs, that keeps us from being lonely and brings us together from different worlds to work on something that’s benefitting people.” – Ronny Hernandez Mora

Round 2 – The Present: What does rOpenSci do well?

The second round invited reflection on what rOpenSci does well today and what participants want to get more involved in.

The most common answer was creating a welcoming space. Participants then commented on the different places where they experience the friendly and welcoming community:

• Co-working came up again here, noted as something other communities have adopted from rOpenSci.

• The Slack workspace and the atmosphere it creates: a common space for learning, getting advice, sharing ideas, and feeling part of something larger.

• Good documentation for package submission, which lowers the barrier to participating.

• Some participants expressed interest in getting more involved in the Champions Program, along with honest questions about whether the entry points are accessible enough.

• The sense that rOpenSci’s quality in different activities comes from its community participation.

You can jump into a role and the community will support you. If there is anything you want to do, start doing it and you will get the support you need.” – Adam Sparks

Round 3 – The Future: What comes next?

The third round asked what participants would like to see rOpenSci do next, and what ideas from other communities could be worth adopting or partnering on.

In-person and hybrid gatherings were a recurring theme: local meetups, local chapters, satellite events at larger R conferences like useR!, and the return of unconferences, potentially co-organised with partner communities.

Other ideas from across the rooms:

• First-time contributor hackathons, possibly with a GenAI focus

• A 24-hour unconference modeled on formats used by other communities like SORTEE

• Career panels showing the range of roles in research software engineering and data-centric work

• An rOpenSci-led package for R package development, like usethis, that supports GitLab, Codeberg, and other platforms beyond GitHub

• Deeper coalitions with open science organizations

• Connecting package authors with graphic designers to help them with their image and visual identity

• Exploring Zarr in R as an area where rOpenSci could help shape the ecosystem

• Better serving very beginners, the people who are still asking “what is R, anyway?”

• Merchandise and swag, we all want our rOpenSci hex sticker

One thread in the conversation pointed to a real tension: rOpenSci has spent 15 years building carefully and steadily. GenAI is moving much faster. How do you maintain that care and those foundations when the pace of change accelerates? Even when we didn’t have any easy answers, this question is the right one to be discussing with the community.

Closing thoughts

Reading through the notes from both events, a few things are consistent.

People talk about rOpenSci not just as a set of tools or infrastructure, but as a place where they found colleagues, had their work taken seriously, and felt like valued contributors.

The software peer-review process, the dev guide, the co-working sessions, the Champions Program, these work because the community behind them is invested in supporting them and find them useful.

At 15 years, the ideas that surfaced – local chapters, hackathons, coalitions, better onboarding for newcomers – suggest a community that still wants to try new things.

It works because everyone is welcome.” – several attendees at the events

Thank you to everyone who joined and shared their thoughts.

These events were part of rOpenSci’s 15th anniversary celebration series. Keep an eye on rOpenSci’s events calendar, our blog and our newsletter for more content around our anniversary.

On 22 July 2026 at 15:05, I made the first commit to a new repository. It held an empty R package and one Excel workbook with field data from a household waste study in Kampala, Uganda. Three minutes later I sent the first prompt to Claude Code:

I am building a new data package and received a dataset in @data-raw/ to process for clean up. […] strictly start with cleaning the data, writing a dictionary. […] make a plan first.

On 24 July at 00:18, the repository held version 1.0.0 of solidwastekampala with a Zenodo DOI: a tidy and validated dataset of 103 households, a data dictionary for 37 variables, documentation, a website, machine-readable metadata, a citation with the ORCID iDs of the authors, and automated checks on five platforms.

My own time on it was no more than about three hours, spread over four sessions. The repository records every step in its commits, issues and pull requests, and in an archive of every prompt I sent.

In version 0.39.0 of the ahead package, the conformalize() function has been updated to include a new parameter called split_prob. This parameter allows users to specify the proportion of the data that should be used for calibration when constructing conformal prediction intervals.

0 – Install packages

install.packages("pak")

pak::pak("Techtonique/ahead")

library(ahead)

1 – USAccDeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = USAccDeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

2 - AirPassengers

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = forecast::ets, y = AirPassengers,
                           h = 20L, split_prob = 0.7)
plot(obj)

3 - fdeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = fdeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

4 - UKDriverDeaths

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = UKDriverDeaths,
                           h = 20L, split_prob = 0.7)
plot(obj)

5 - WWWusage

par(mfrow=c(2, 2))
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.4)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.3)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.6)
plot(obj)
obj <- ahead::conformalize(FUN = ahead::dynrmf, y = WWWusage,
                           h = 20L, split_prob = 0.5)
plot(obj)

Thanks to the awesomeness of webR and Shinylive, the tool I made in back in 2014 to estimate your blood alcohol concentration (BAC) now runs 100% in the browser, no server needed!

Try the browser version of drinkR here!

library(knitr)
library(ggplot2)
library(dplyr)
library(tidyr)
library(broom)

set.seed(123)
n <- 100
X1 <- rnorm(n, mean = 0, sd = 1)
X2 <- rnorm(n, mean = 0, sd = 1)
Z <- 5 * X2 + rnorm(n, mean = 0, sd = 1)
Y <- 3 * X1 + 2 * X2 + rnorm(n, mean = 0, sd = 1)

kable(cor(cbind(X1, X2, Z)), caption="Correlation of the predictors")

fit1 <- lm(Y ~ X1 + X2)
fit2 <- lm(Y ~ X1 + X2 + Z)
coef_results <- rbind(
  tidy(fit1) %>% mutate(model = "fit1"),
  tidy(fit2) %>% mutate(model = "fit2")
)
# move model to the first column
coef_results <- coef_results %>% select(model, everything())
# now R-squared and adjusted R-squared
rsq_results <- rbind(
  glance(fit1) %>% mutate(model = "fit1"),
  glance(fit2) %>% mutate(model = "fit2")
) %>%
  select(model, everything())
kable(coef_results, caption = "The coefficients and statistics")

kable(rsq_results, caption = "The goodness of fit")

simulate_correlated_data<-function(rho, n=100, seed=NULL) {
  if(is.null(seed)) {
    seed < Sys.time()
  }
  set.seed(seed)
  x <- rnorm(n)
  z <- rnorm(n)
  y <- rho*x + sqrt(1-rho^2)*z
  return(cbind(x,y))
}

n<-100
rhos<-c(0.1, 0.3, 0.5, 0.8, 0.95)
simulated_data <- lapply(rhos, simulate_correlated_data, n=n)
# let's calculate the correlation coefficients for the simulated data
simulated_rhos<-sapply(simulated_data, function(d) cor(d[,"x"],d[, "y"]))
kable(cbind(expected_rho=rhos, actual_rho=simulated_rhos))

 fits <- lapply(simulated_data, function(d) {
          x1<-d[,"x"]
          x2<-d[, "y"]
          y <- 10*x1 + 3*x2 + rnorm(n)
          fit<-lm(y~x1 + x2)
          return(fit)
          } 
          )
coef_results<-lapply(seq_along(rhos),
                     function(i) {
                       tidy(fits[[i]]) %>%
                         mutate(rho = rhos[i]) %>%
                         select(rho, everything())
                     }
                     )
coef_results <- do.call(rbind, coef_results)
kable(coef_results, caption = "Regression statistics with different correlation strengths among predictors")

coef_results <- coef_results %>% 
    mutate(rho = factor(rho)) %>%
    pivot_longer(
      cols = estimate:p.value,
      names_to = "metric",
      values_to = "value"
    )
plt <- ggplot(coef_results, aes(x=rho, y=value)) +
      geom_point(aes(color=term)) +
      facet_wrap(~ metric, ncol=2, scales="free_y") +
      theme(legend.position = "top")
plt

If you plan to learn R or Python in the new year, book before the end of November. Our spring 2027 data science training courses run from January to June, and early bird pricing on the first one ends on 29 November.

There are 22 courses on the schedule. All of them run live online, and the trainers are data scientists and engineers who spend the rest of their week on client projects. Each course runs over two afternoons, 13:30 to 17:00 UK time. The Bayesian inference course is the exception: it runs over four.

Three routes through our spring 2027 training courses

The courses are scheduled so you can take them in order. There are three obvious routes, each ending with a set of skills you can use at work the following week.

Learning R from scratch

Start with Introduction to R on 11–12 January. Data Wrangling in the Tidyverse follows on 25–26 January, then Data Visualisation with ggplot2 on 1–2 February and Programming with R on 8–9 February. That gets you from your first script to writing your own functions in a month.

Why take this route: Clean, chart and automate your data in code. Monthly reporting becomes a script you rerun in minutes, with fewer manual errors.

January too soon? The same four courses run again between 12 April and 11 May. If you are still deciding whether R is worth the time, we wrote about why learning R is a good career move.

Learning Python

Introduction to Python runs on 15–16 February. Programming with Python follows on 22–23 February, and Data Visualisation with Python on 1–2 March. The same three run again in May and June.

Why take this route: Analyse data in pandas, automate repeat tasks and produce report-ready charts, in the same language as your engineering team. It is also the groundwork for machine learning in Python.

Statistics and machine learning in R

This route is for people who already use R. Statistical Modelling with R on 8–9 March covers hypothesis testing, regression, clustering and principal components analysis. Machine Learning with Tidymodels on 15–16 March and Advanced Machine Learning with Tidymodels on 22–23 March take you through pre-processing, model fitting and tuning.

Why take this route: Build, validate and tune predictive models in tidymodels, and explain which variables drive them. Your team moves from reporting what happened to forecasting what happens next, with evidence for how far to trust it.

Two courses that run once this spring

Introduction to Bayesian Inference using RStan runs over four afternoons, 18–21 January. You write Stan programs for a range of statistical models and learn to read the MCMC diagnostics in R. It is the only date for it in the first half of 2027.

Reporting with Quarto runs on 21–22 June. You build reports, slides and documents straight from R or Python code, so the numbers in the report change when the data does. It pairs well with the ideas in our post on reproducible analytical pipelines.

Early bird deadlines

Early bird pricing ends about six weeks before each course. Booking closes a week before the start date. The first three deadlines are:

• 29 November for Introduction to R
• 6 December for Bayesian inference
• 13 December for Data Wrangling in the Tidyverse

Dates for the second half of 2027 will follow.

Using this year’s training budget

If your team has training budget left this December, booking a January course now spends it before the year ends. It also gets you the early bird price, because every January deadline falls before Christmas.

What you get on every course

• PDF notes and scripts to keep
• Live demonstrations and hands-on exercises
• Time with the trainer to ask questions about your own work

We have run more than 1,000 courses, for organisations including NHS Scotland, Shell, Wessex Water and the Royal Statistical Society.

Training a whole team

Every course in our full catalogue can also run privately for one team, using your data and your tools. Group bookings and returning clients get a discount. Email training@jumpingrivers.com with the size of the team and what they need to learn.

See the full public schedule and book a place.

For updates and revisions to this article, see the original post

The Model Context Protocol (MCP) is quietly becoming one of the most useful ideas in AI tooling. Here, we describe what MCP is, why it matters, and how to build an MCP server yourself.

MCP is a standard way to link a Large Language Model (LLM) to the outside world. An MCP server is a small piece of software that exposes some data or capability to an AI model. Everyday examples include connecting your AI agent to your calendar or your file system, or allowing it to search the web. These are things an LLM cannot do (or would struggle to do) by itself. On its own, a model only knows what it learned during training, plus whatever you provide it, and it has no way to check, fetch or act on anything beyond that. Even then, it can only give out text responses, which you have to be the one to act on. An MCP server can allow an LLM to carry out actions itself. Less talking, more doing.

Why not just use an API?

There are many ways to connect external capabilities to an AI model, usually via an API (Application Programming Interface), and this was the industry standard. However, APIs are all configured slightly differently, and once you start to scale up the number of connections you need between multiple pieces of software, you hit a combinatorial problem: the number of custom connections a product needs multiplies by the number of tools that use it.

Anthropic introduced MCP in late 2024 as an open standard for connecting tools and data to LLM applications. It works like a USB port on a computer: any MCP server can be used with any agent/model, whether you are using Claude, Gemini, ChatGPT or one you’ve made yourself.

Building your own MCP server

Full instructions on how to build an MCP server are available in the MCP documentation, but here is a very quick overview using Claude Code and the Sonnet 5 model.

The good news is that a working MCP server can be as short as a few lines of code. Here’s the shape of it, step by step, using Python.

Step 1: Install the dependencies

The MCP project maintains official SDKs for Python, TypeScript and other languages. The Python ‘mcp’ package embodies this SDK. However, here we are going to use the ‘fastmcp’ package. This is built upon ‘mcp’, but has some benefits (in particular, ‘fastmcp’-based tools can readily be deployed on Posit Connect, see later).

mkdir -p ~/projects/
cd ~/projects
uv init --python=3.15 dice-roll
# Initialized project `dice-roll`
# at `/home/<user>/projects/dice-roll`
cd dice-roll
uv add fastmcp

This gives us the following project structure, which is typical for a ‘uv’ project. We will define our MCP server in ./src/dice_roll/__init__.py, although if the project got much larger, we would split it into separate modules or sub-packages.

tree -L 2
.
├── pyproject.toml
├── README.md
├── src
│   └── dice_roll
└── uv.lock

Step 2: Define a tool

An MCP server exposes “tools”: functions the model can call. Here is a mini example: a server with a single tool that rolls a dice.

# src/dice_roll/__init__.py
from typing import Annotated
import random

from fastmcp import FastMCP
from pydantic import Field

# Give your server a name: this is what shows up when a model connects
mcp = FastMCP("dice-pick")

@mcp.tool
def roll_dice(
 sides: Annotated[int, Field(description="Number of sides on the dice", ge=2)] = 6,
) -> int:
 """Roll a single fair dice with the given number of sides and return the number rolled."""
 return random.randint(1, sides)

def main() -> None:
 mcp.run()

The @mcp.tool docstring (“Roll a single …”) is important. This becomes the description of the tool that the LLM reads. It is like the blurb of a book: the LLM will look at the description and from it decide whether it wishes to use the tool or not. If this description isn’t clear, it may lead to the tool not being used.

The type annotations for both the parameter sides and the value returned by roll_dice() also provide useful metadata for the LLM. Our sides parameter has a constraint that ensures at least 2 sides are present on the dice (ge=2; this is enforced within Python by Pydantic). When describing the tools provided by this MCP server, an LLM will be able to explain that the minimum number of sides is 2 (and more advanced models may not call the tool otherwise).

Step 3: Connect it

Save the code in the project, and we will test whether an MCP-compatible agent can use it. There are a few ways to register our locally-defined MCP server with Claude Code.

Since we have used the fastmcp framework, there is a command line tool that installs the MCP server for use from the project directory:

# From the ~/projects/dice-roll/ directory
uv run \
 fastmcp install claude-code src/dice_roll/__init__.py:mcp --with-editable .

The uv run part is necessary to ensure the correct Python environment is used. The fastmcp command needs to find the mcp object that we have defined in our source code. If we open claude and ask “What dice-pick tools do you know about?”, it replies:

The dice-pick MCP server gives me one tool:
roll_dice: rolls one fair dice and returns the number it lands on.
┌───────────┬─────────┬─────────┬───────────────────┐
│ Parameter | Type │ Default │ Constraint │
├───────────┼─────────┼─────────┼───────────────────┤
│ sides │ integer │ 6 │ must be 2 or more │
└───────────┴─────────┴─────────┴───────────────────┘
So I can roll a d6, d20, d100 or any other single dice. It can't roll several dice in one call, add
modifiers or pick from a list. For something like 3d6, I'd call it three times.

Then we can use the tool:

❯ Roll a dice
Called dice-pick
I rolled a 6 on a standard six-sided dice.
...
❯ Roll a 20-side dice
Called dice-pick
I rolled a 20 on the d20.
That's the highest number on both rolls so far. [... blah blah blah ...]
...
❯ Roll a d12
Called dice-pick
I rolled a 4 on the d12, so the run of top scores is over.

We have installed our locally in-development MCP server on our machine, but it’s only usable from the ~/projects/dice_roll directory. This is fine while you are testing, but you may want to use your MCP server in other locations once it works well.

Our project is uv-based. If we tell it which module and function should be called to start the project running, we can make the project more easily usable from other locations. But, the entry-point function used to start our MCP server running is dice_roll.main(), which uv has already added as the entry-point to the project in the pyproject.toml:

[project.scripts]
dice-roll = "dice_roll:main"

On running uv run python or uv sync, a dice-roll entry-point will be added to the virtual environment for the project.

For Claude Code, we can then register our MCP server with a single command:

claude mcp add --scope user dice-pick \
 -- /abs/path/to/projects/dice-roll/.venv/bin/dice-roll

Now, you can ask claude to “roll a 20-sided dice” from any directory on your computer and it will call your tool directly rather than guessing at an answer. This is possible because claude mcp add ... appends to the mcpServers entry of a file (~/.claude.json).

A few tips before you build something bigger:

• Keep tool descriptions (the docstrings for your @mcp.tool-decorated functions) short and specific. The model relies on the docstring to decide when to use a tool, so vague descriptions lead to it being ignored or misused.
• Think about what should require confirmation. Read-only tools (like looking something up) are usually safe to let a model call freely; tools that change or delete something are worth gating behind an explicit user confirmation step.

You can remove the MCP server from your claude setup by using

claude mcp remove dice-pick

Or just ask Claude to remove it!

Deploying to Posit Connect

A server running on your own machine is fine for testing, but it is only useful to you. MCP servers can be deployed to a range of different cloud providers. Posit Connect, the platform we normally reach for when publishing Shiny apps, Plumber APIs and Quarto reports, added first-class support for hosting MCP servers in 2026, complete with a dedicated icon and its own tab in the content list.

Posit Connect does expect you to write your Python-based MCP servers using the fastmcp framework that we used above. Connect also needs to import a runnable application object that can handle HTTP-based requests. When deployed, ‘fastmcp’-based MCP servers can expose an ASGI app object that handles these requests. Locally run servers use ‘stdio’ routing, but we only have to modify our current source code slightly to deploy to Posit Connect.

# imports ...

mcp = ...

def roll_dice():
 ...

app = mcp.http_app() # For deployment on Connect

def main():
 mcp.run() # For running using stdio locally

For deployment of the MCP server to Posit Connect, we also need to create a requirements.txt file. These are used to specify the run-time (and sometimes development) dependencies of a project. Since we used uv to manage the project, this is easy to create:

uv export --no-dev --no-hashes -o requirements.txt

The requirements.txt file contains a line -e . that dynamically includes any source code from src/. The Python environment that is built for your MCP server depends upon requirements.txt and is rebuilt when that file is updated.

We use the rsconnect-python package (and the rsconnect command line tool that it installs) to interact with Posit Connect. You may first have to register with your Connect system using rsconnect add ... or rsconnect login ....

Then you can deploy with the rsconnect tool, using the FastAPI deployment mode and pointing the entrypoint at the app object that we have just added to the source code. This is done using a --entrypoint argument of the form module:object_name. Here this is dice_roll:app, since our app object is in the dice_roll module:

rsconnect deploy fastapi -n <nickname-of-posit-server> --entrypoint dice_roll:app ./

The Posit Publisher extension for VS Code and Positron works too, if you’d rather deploy from your editor. Posit has some good documentation on publishing MCP servers.

On first deployment, Connect scans requirements.txt, recognises fastmcp, and automatically sets the content type to MCP. Once deployed, it keeps at least one process warm so there is no cold-start delay, and gives it its own icon in the content list, with no extra configuration needed.

Using your deployed MCP server

To use an MCP server that is deployed in the cloud, the process is very similar to how we added our locally-defined MCP server. On Posit Connect (or wherever you have published), you will have a URL for the MCP server and you can configure claude to use this URL as follows:

claude mcp add --transport http <mcp-name> https://<connect-server>/content/<content-hash>/mcp

Here we used dice-pick as the name of the MCP server, consistent with the name in the source code. When you next open claude you will be prompted to authenticate before you can use the tool.

Want to see this in action? Watch this walkthrough, recorded live at our AI in Production 2026 conference, of an MCP server being built and used from scratch.

Key takeaways

• An MCP server is the standard way to give an AI model access to things it can’t do alone: checking a calendar, reading a file, searching the web, or looking up live data.
• It works like a USB port: build the connection once, and any compatible model can use it.
• A working MCP server can be as short as a few lines of code. The model calls your tools directly instead of guessing at an answer.
• Keep tool descriptions short and specific, and gate anything destructive behind explicit confirmation.
• Posit Connect can host your MCP server directly, so your whole team can use it rather than just your own machine.

For updates and revisions to this article, see the original post

One of the fundamental aims of field experiments is to explorare the relationships between several experimental traits, which have been recorded/measured altogheter for all the experimental units. For example, think about a genotype experiment with 27 winter wheat genotypes in a Randomised Complete Block Design with three replicates. If we aim at assessing yield and yield quality, we have to measure, in each plot, several traits, such as grain yield and the weight of thousand kernels (TKW). Apart from separately analysing these two variables, we might be interested in describing their reciprocal relationships, by using the Pearson’s correlation coefficient. Let’s do this, by using the dataset ‘WheatQuality’, which is available in the ‘statforbiology’ package; we have 81 plots and just as many couples of measures in all.

library(Hmisc)
library(dplyr)
library(statforbiology)
dataset <- getAgroData("WheatQuality")
head(dataset)
##     Genotype Block Height  TKW Whectol Yield
## 1 arcobaleno     1     90 44.5    83.2 64.40
## 2 arcobaleno     2     90 42.8    82.2 60.58
## 3 arcobaleno     3     88 42.7    83.1 59.42
## 4       baio     1     80 40.6    81.8 51.93
## 5       baio     2     75 42.7    81.3 51.34
## 6       baio     3     76 41.1    81.1 47.78

Now, we consider the whole lot of measures and calculate the correlation coefficient between yield and TKW; we also use a t-test for the hypothesis that the correlation coefficient is equal to 0. This is the result:

ra <- with(dataset, rcorr(Yield, TKW) )
# 
# Correlation coefficient
ra$r[1,2] 
## [1] 0.540957
#
# P-level
ra$P[1,2] 
## [1] 1.850931e-07

The conclusion is: the correlation coefficient is positive and significantly different from 0, which supports the idea that plots with a high value on yield tend to have a high value on TKW, as well.

Such an approach is very often seen also in high ranking Journals, but: is it correct? Can we trust the conclusions? The answer is no, not completely. And here is why. Please, note that in this post I am following the point made by Bland and Altman (1994), by using an example that is more closely related to plant/crop science.

# Correlations within genotypes
wgCors <- by(dataset, list(dataset$Genotype), function(x) cor(x$Yield, x$TKW))
wgCors
## : arcobaleno
## [1] 0.9847228
## ------------------------------------------------------------ 
## : baio
## [1] 0.1611952
## ------------------------------------------------------------ 
## : claudio
## [1] -0.9993872
## ------------------------------------------------------------ 
## : colorado
## [1] 0.9837293
## ------------------------------------------------------------ 
## : colosseo
## [1] 0.4564855
## ------------------------------------------------------------ 
## : creso
## [1] -0.5910193
## ------------------------------------------------------------ 
## : duilio
## [1] -0.988233
## ------------------------------------------------------------ 
## : gianni
## [1] -0.7603802
## ------------------------------------------------------------ 
## : giotto
## [1] 0.9520211
## ------------------------------------------------------------ 
## : grazia
## [1] 0.4980828
## ------------------------------------------------------------ 
## : iride
## [1] 0.7563338
## ------------------------------------------------------------ 
## : meridano
## [1] 0.1174342
## ------------------------------------------------------------ 
## : neodur
## [1] 0.4805871
## ------------------------------------------------------------ 
## : orobel
## [1] 0.8907754
## ------------------------------------------------------------ 
## : pietrafitta
## [1] -0.9633891
## ------------------------------------------------------------ 
## : portobello
## [1] 0.9210135
## ------------------------------------------------------------ 
## : portofino
## [1] -0.9900764
## ------------------------------------------------------------ 
## : portorico
## [1] 0.1394211
## ------------------------------------------------------------ 
## : preco
## [1] 0.9007067
## ------------------------------------------------------------ 
## : quadrato
## [1] -0.5840238
## ------------------------------------------------------------ 
## : sancarlo
## [1] -0.646067
## ------------------------------------------------------------ 
## : simeto
## [1] -0.4051779
## ------------------------------------------------------------ 
## : solex
## [1] -0.6066363
## ------------------------------------------------------------ 
## : terrabianca
## [1] -0.4076416
## ------------------------------------------------------------ 
## : verdi
## [1] 0.5801404
## ------------------------------------------------------------ 
## : vesuvio
## [1] -0.7797493
## ------------------------------------------------------------ 
## : vitromax
## [1] -0.8056514

# Correlations within blocks
wbCors <- by(dataset, list(dataset$Block), function(x) cor(x$Yield, x$TKW))
wbCors
## : 1
## [1] 0.5998137
## ------------------------------------------------------------ 
## : 2
## [1] 0.539999
## ------------------------------------------------------------ 
## : 3
## [1] 0.5370398

mod1 <- lm(Yield ~ factor(Block) + Genotype, data = dataset)
mod2 <- lm(TKW ~ factor(Block) + Genotype, data = dataset)
wCor <- rcorr(residuals(mod1), residuals(mod2))
wCor$r
##            x          y
## x 1.00000000 0.03133109
## y 0.03133109 1.00000000
wCor$P
##           x         y
## x        NA 0.7812693
## y 0.7812693        NA

means <- aggregate(dataset[,c(4,6)], list(dataset$Genotype), mean)
rgPear <- rcorr( as.matrix(means[,2:3]) )
rgPear$r
##             TKW     Yield
## TKW   1.0000000 0.5855966
## Yield 0.5855966 1.0000000
rgPear$P
##              TKW      Yield
## TKW           NA 0.00133149
## Yield 0.00133149         NA

means <- aggregate(dataset[,c(4,6)], list(dataset$Block), mean)
rbPear <- cor( as.matrix(means[,2:3]) )
rbPear
##               TKW       Yield
## TKW    1.00000000 -0.08812544
## Yield -0.08812544  1.00000000